function [B,G] = inflow_bc(V,T,TE,ET,d,dof_map,cr,bdr,desc,Re,caseNum)

% function [B,G] = inflow_bc(V,T,TE,ET,d,dof_map,cr,bdr,desc,Re,caseNum)

% this is the inflow boundary condition

% the common boundary can be viewed as inflow boundary conditions

%

nbdr = length(bdr);

Indx2 = zeros(nbdr*(2*d+1),1);

G = zeros(nbdr*(2*d+1),1);

row = 1; pos = 1;

G0 = G(1);

% begin to treat boundary condition one by one bound edge

for k = 1:nbdr

    eg = bdr(k);  tri = ET(eg,1);  

    v1_loc = find(TE(tri,:)==eg);

    v2_loc = mod(v1_loc,3)+1;

    v3_loc = mod(v2_loc,3)+1;

    V1 = V(T(tri,v1_loc),:); 

    V2 = V(T(tri,v2_loc),:);  

    V3 = V(T(tri,v3_loc),:);  

    row_idx = row:row + 2*d;

    if v1_loc==2;

        v1_loc = -2;

    end    

    % in tcord, the vertices order must use physical order

    line1 = cr_indices(0,d-1,v1_loc,cr);

    A = [1,1,1;V(T(tri,1),1),V(T(tri,2),1),V(T(tri,3),1);...

        V(T(tri,1),2),V(T(tri,2),2),V(T(tri,3),2)];

    x = A\[0;1;0];

    mat_x = desc_mat(d,x(1),x(2),x(3),desc);

    y = A\[0;0;1];

    mat_y = desc_mat(d,y(1),y(2),y(3),desc);

    Dx = mat_x(line1,:);

    Dy = mat_y(line1,:);

    [i,j,s] = find([Dx;Dy]);

    L = length(i);

    Indx2(pos:pos+L-1) = dof_map(j,tri);

    pos = pos + L;

    %%%%%%%%%%%%%%%%%%%%%%% not finished yet%%%%%%%%%%%%%%%%%%%%

    line0 = cr_indices(0,d,v1_loc,cr);

    line1 = cr_indices(1,d,v1_loc,cr);

    [c1,c2] = navstk_first_layer(G0,V2,V3,V1,d,Re,caseNum); 

    G(row_idx) = [c1;c2];  

    Indx2(row_idx) = dof_map([line0;line1],tri);

    row = row + 2*d + 1;

    G0 = c1(d+1);

end

ONE = ones(row-1,1);

Indx1 = (1:row-1)';

dim = max(max(dof_map))-min(min(dof_map)) + 1;

B = sparse(Indx1,Indx2(1:row-1),ONE,row-1,dim);

G = G(1:row-1);

end



function [c1,c2] = navstk_first_layer(G0,V1,V2,V3,d,Re,caseNum) 

d = d-1;

I = (0:d)';

J = d - I;

X = (J*V1(1) + I*V2(1))/d;

Y = (J*V1(2) + I*V2(2))/d;

[val,G1,G2] = boundary_function(X,Y,Re,caseNum);

b3 = -G2*(V3(1)-V1(1)) + G1*(V3(2)-V1(2));

b2 = -G2*(V2(1)-V1(1)) + G1*(V2(2)-V1(2));

m = (d+1);

IM = diag(I)*ones(m,m);

JM = diag(J)*ones(m,m);

Mat = (IM/d).^(IM').*(JM/d).^(JM');

IF = gamma(I+1);

JF = gamma(J+1);

A = factorial(d)*ones(m,m)*diag(1./(IF.*JF));

Mat = A.*Mat;

D2 = 1/(d+1)*(Mat\b2); % D2 is the value of c^{(1)}

B = diag(ones(d+1,1))-diag(ones(d,1),-1); % this is the mat of D_d(bary(\xi)), 

                                % actually bary(\xi) = (-1,1,0) here

D2(1) = D2(1)+G0;

c1 = B\D2; % get the d B-coefficient except for the first one. this may be not accurate , I add myself

%%%%%%%%%%%%%%%%%% notice here%%%%%%%%%%%%%%%%%%%%%%

c1 = [G0;c1]; % add the first one

c2 = 1/(d+1)*(Mat\b3) + c1(1:(d+1)); % comput the second line on boundary, here n3=1,n2=0,n1=-1

end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%



function line_dof = cr_indices(r,d,idx,pattern)

% to get the r's line of dofs of degree d parallel to edge i

% it is specially useful when treating smoothness condition and boundary

% condition! 

% the output parameter mat is the matrix P in my notes.

switch idx

    case 1

        line_dof = (((d+1-r)*(d-r)/2+1):((d+2-r)*(d+1-r)/2))';

    case -1

        line_dof = (((d+2-r)*(d+1-r)/2):-1:((d+1-r)*(d-r)/2+1))';

    case 2

        line_dof = pattern(r+1,1:d-r+1)';

    case -2

        line_dof = pattern(r+1,d-r+1:-1:1)';

    case 3

        line_dof = pattern(1:d-r+1,r+1);

    case -3

        line_dof = pattern(d-r+1:-1:1,r+1);

    otherwise

        line_dof = [];

end

%then put the corresponding value to a matrix;

% i = (1:d-r+1)';j = line_dof; s = ones(d-r+1,1);

% mat = sparse(i,j,s,d-r+1,(d+1)*(d+2)/2);

end



%cccc

function desc = desc_mat(d,lam1,lam2,lam3,desc_pattern)

% get the corresponding matrix for evaluating algorithms

% from degree d to d-1

m_rows = d*(d+1)/2;

m_cols = (d+1)*(d+2)/2;

I = (1:m_rows)';

Id = ones(m_rows,1);

desc = sparse(I,desc_pattern(I,1),lam1*Id,m_rows,m_cols) + sparse(I,desc_pattern(I,2),lam2*Id,m_rows,m_cols) + sparse(I,desc_pattern(I,3),lam3*Id,m_rows,m_cols);

end